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Abstract 

We consider four turbulent models to simulate the boundary mixing layer of the 
ocean. We show the existence of solutions to these models in the steady-state case 
then we study the mathematical stability of these solutions. 
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1 Introduction 

The presence of an homogeneous layer near the surface of the ocean has been observed since 
a long time. The so called "mixed layer" presents almost constant profiles of temperature 
and salinity (or equivalently the density). The bottom of the mixed layer corresponds 
either to the top of the thermocline, zone of large gradients of temperature, or to the top 
of the zone where haline stratification is observed [5]. Some attempts to describe this 
phenomenon can be found for example in Defant [3] or Lewandowski [5j. The effect of the 
wind-stress acting on the sea-surface was then considered to be the main forcing of this 
boundary layer. Observations in situ were completed by laboratory experiments [2] and 
more recently by numerical modelizations of the mixed layer. 

In this note, we consider four turbulent models to describe this homogeneous boundary 
layer. The first one is the Pacanowski-Philander model, and two of these models are new 
models. They aim to compute the velocity and the water density of a water column, 
are one space dimensional and the eddy viscosities depend on the Richardson number. 
For those model, we show the existence of a steady-state solution and we analyse the 
mathematical linear stability of these steady state solution, showing that only one of these 
model, the one we introduce in this note (model labelized as R — 2 — 2 — 4 below), has a 
unique staedy state solution with a large range of stabilty. Moreover, in [I] we have used 
these models to simulate the warm pool at the equator. Numerical results confirm that 
R — 2 — 2 — 4 is the most accurate parametrization. 



2 The equations 

We denote by (u, v) the horizontal water velocity and p its density. Since the numerical 
simulation performed in [1] concerns the equator zone, we do not take the Coriolis force 
into account. The closure equations are: 
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u = uq, v = Vq, p = po at initial time t = 0. 



In system (|2.ip , the coefficients ui and ^2 are the vertical eddy viscosity and diffusivity 
coefficients and will be expressed as functions of the Richardson number R defined as 
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where g is the gravitational acceleration and po a reference density (po ~ 1025 kg.m 3 ). 

The constant h denotes the thickness of the studied layer that must contain the mixing 
layer. Therefore the circulation for z < —h, under the boundary layer, is supposed to be 
known, either by observations or by a deep circulation numerical model. This justifies 
the choice of Dirichlet boundary conditions at z = —h, Ub, Vb and pb being the values 
of horizontal velocity and density in the layer located below the mixed layer. The air-sea 
interactions are represented by the fluxes at the sea-surface : V x and V y are respectively the 
forcing exerced by the zonal wind-stress and the meridional wind-stress and Q represents 
the thermodynamical fluxes, heating or cooling, precipitations or evaporation. We have 
V x = Co \u a \ 2 and V y = Cd \v a \ 2 , where U a = (u a ,v a ) is the air velocity and Cd a friction 
coefficient. 

We study hereafter four different formulations for the eddy coefficients V{ = fi (R) , 
labeled as "R - 2 - i" and/or "R-2-i-f . In all models, /1 (R) = a x + 
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Formulation 12 — 2 — 1 — 3 corresponds to the modelization of the vertical mixing proposed 
by Pacanowski and Philander [7j. The coefficients ai,/?iand a 2 have the following values: 



a\ = 1.10 -4 , Pi = 1.10 -2 , «2 = 1.10 _5 (units:m 2 s _1 ).This formulation has been used 
in the OPA code developed in Paris 6 University [B] with coefficients a\ = 1.10 -6 , j3\ = 
1.10 -2 , «2 = 1-10 -7 (units: m 2 s _1 ) .The selection criterion for the coefficients appearing 
in these formulas was the best agreement of numerical results with observations carried 
out in different tropical areas. Formulation R — 2 — 3 has been proposed by Gent [I]. 
Formulations R — 2 — 2 — 4 and R — 2 — 2 are new as far as we know. Notice that models 
R — 2 — 1 — 3 and R — 2 — 3 are no more physically valid respectively for R G (—3.13, —0.2) 
and R £ (—2.25, —0.1) since the coefficient v<i becomes negative. 

2.1 Steady-state solutions 

Steady-state solutions to system (|2.ip satisfy 
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Theorem 2.1 System ([23$ has at leat one smooth solution on [0, — h] for each model in 
(\2.2§ . In case of R — 2 — 2 — 4 £/ie solution is unique. 

Proof. Integrating ()2.3[) with respect to z, yields 
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which is a fixed point equation for R. Any solution R to equation (12, 5p yields a Richardson 

number i? e corres ponding to the fluxes V x ,V y and Q and not on z as i^i and z/2 are 
independent on the depth variable z as well asq the turbulent viscosities. The Richardson 
number R e being known, steady-state profiles for velocity and density are obtained by 
integrating (|2.4p with respect to z, taking into account the boundary conditions at z = —h: 

u* (z) =u b + V f a (z + h), v*{z)=v b + V f \ ' (* + /*), 
(2.6) P°j£ ( R > P°J l y R > 

pe{z) = pb+ MR T ) {z + h) - 

It remains to analyse the existence of solutions of equation (j2.5|) . These solutions can 
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be interpreted as the intersection of the curves k (R) = w \ and h (R) = CR with 

P 2 a(V x 2 + V v 2 ) 

C = — . The existence and the number of solutions are controlled by the 
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constant C and then by the parameter — — , V 2 = V 2 + V 2 , depending only on the surface 

fluxes. The graph of function k and h for Q < and Q > is plotted on Figures 1 and 2 
below when fi and fi in case of R-2-2-4 and R — 2 — 2. 




Formulation R-2-2-4 Formulation R-2-2 

Figure A Figure B 

The qualitative behaviour obtained with formulation R-2-3 and R-2-1-3 is the same as 
R-2-2. The intersection of k(R) and h(R) = CR consists in one point for Q < and 
several points for Q > 0. The number of points depends to the values of surface fluxes. 

The graphs obtained for the R-2-2-4 modelization (Figure A) and its simplified version 
R-2-2 (Figure B) are very different. It is obvious in Figure A that any straight line h (R) = 
CR meets k at only one point for Q > and Q < 0. Therefore it exists one unique 
equilibrium Richardson number R e whatever the values of the surface fluxes V x , V y and 
Q. In the case of the other models, we get several solutions. The proof is finished. Notice 
that in [1] we show that the most accurate model is R — 2 — 2 — 4 from the physical and 
numerical viewpoint. 



2.2 Linear stability of the equilibrium solutions 

In this section we analyse the time evolution of a small perturbation of one of the equilib- 
rium states (u e ,v e ,p e ) described in the previous section. 

At initial time t = we set (uq,vq,pq) = (u e ,v e ,p e ) + (u' ,v' ,p' ) and we denote by 

(u,v,p) = (u e ,v e ,p e ) + [u',v',p') 

the solution of equations (|2,ip at time t where (u e ,v e ,p e ) are solution to the steady-state 
system (|2.3p . and v\ = f± (R e ) and = f2 (R e ) are two positive constants. 

d p du dv 
Introducing the new variables tp = — , 9 = — — and (3 = — , the Richardson number 
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can be expressed as 

Po (0 2 + p 2 ) 
Applying the Taylor formula, we get 
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We set for k = 1, 2 : T = v\ (9, [3, ip) - (0 e , /? e , ^ e ), 9 = ^2 (0, /?, V) - ^2 (# e , /? e , V> e ) and 
*/| = v k {9 e ,p e ,ip e ), 9' = 0-e e , (5' = (3-(5 e , ip' = ip~ip e , 



( duk \ 

\d9j 



du k 



89 



(9 e ,/3 e ,r 



dv k 



(0 e ,f3 e ,iP e 



dip) 



dv k 



dip 



' dv k 
dp J dp 

The equations satisfied by the perturbation (u',v',p') are deduced from equations (|2.ip : 
du' d , „ dv' d 
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We now replace v\ and z/2 by expresions deuced from the Taylor's development and retain 

only the first order terms. The approximated equations for (u',v ,p') then are 
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Equations (|2.8|) can be written 
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Let (Ai , A2, A3) be the eigenvalues of matrix A. Assuming the eigenvalues dinstincts, matrix 
A is equal to P~ 1 DP, where D is diagonal, and such that du = Ai, c?22 = A2 and e?33 = A3. 

dW d 2 W 

Set now W = PV. The vector W verifies the system — D = 0, i.e. 



dt 



dz 2 



(2.10) ^l-xf Wi 



dw 2 d 2 w 2 dw 3 d 2 w 3 
U, —^r. a 2 - = U, — — A3- 



0. 



dt 1 dz 2 dt dz 2 dt dz 2 

Stability of the equilibrium solution (u e , v e , p e ) means that any perturbation (u' , v' , p' ) 
imposed at initial time t = is damped as t — > 00. This is verified if the eigenvalues 
Ai, A2, A3 are such that Re (Ai) > 0, Re (X 2 ) > and Re (A3) > 0. These three conditions 
are equivalence to det (A) > 0, tr (A) > and tr (Adj(A)) > 0. From these conditions, we 
build the graph below (see figure 1), obtained thanks an analytical computation (we skip 
the technical details here): 



The results are summarized in Figure 1. 
The circle zone represents a zone where 
the solution is physically not valid. It 
is the case for the R-2-3 and R-2-1-3 
formulation. The rectangular zone is a 
unstability zone. All formulations have 
a unstability zone. Nevertheles, one ob- 
serves that for each model, mathemat- 
ical stability holds for non negative R. 
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Figure 1: Numerical stability 

3 Conclusion 

All the models have a steady-state solution, unique in the case of R-2-2-4. Each one 
is linearly stable for non negative R, which corresponds to physical stability. All these 
models present a linear unstable zone, located in a region where R is non positive. They all 
presents a linear stability zone for some non positive values of R, situation that can arise in 
real situation, as reported in pQ (physical unstability). All these models have been tested 
in pp. The simulation confirms the existence of stable linear steady-state solutions and 
the ability of these models to describe a boundary mixing layer. However, the numerical 
study in pQ confirms that R — 2 — 2 — 4 yields better numerical results. 
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